!##############################################################################
! MAIN FILE 
!
! This code is the main file to solve the life-cycle problem with 
! couples, singles, and portfolio choice
!
! Author: Annika Bacher (annika.bacher@bi.no)
!          
! This version: May 2024
!
!##############################################################################

include "globals.f90"
include "retirement.f90"
include "simulate_economy.f90"

program main

        use globals
        use toolbox_own
        use simulate_economy
        use retirement

        ! declare private variables in main file
        
        implicit none
        
        real*8  :: maxbndm, minbndm, maxbndf, minbndf, maxbndc, minbndc,  output_f(5), output_m(5), output_c(5)
        
        integer :: inx, inx2(1)
        
        integer :: t, p, z, zp, tp, ll

        ! specify parameters 
        
          betaf  = 0.9168d0    
          betam  = betaf  
          betac  = betaf 	                          
          gammaf = 5.0575d0 	      
          gammac = gammaf				   
          gammam = gammaf  
          
        ! age-dependent participation costs
         
		  SF_young = 0.0576d0 
          SF_old = SF_young + 0.06d0   
          SMFs(1:16) = SF_young 
          
          numb = JJ-17
          step = (SF_old - SF_young) / numb
  
		  do i = JJ-numb, JJ
		  SMFs(i) = SF_young + (i - JJ+numb) * step
		  end do 
       
		  SMFc = SMFs 	
          
          p_rec  = 0.146d0	! i.i.d probability of recession 	    
          
        ! return on safe asset

          r = 0.020d0   
          
        ! initialize variables
            call initialize() 
            
        ! continuation values during retirement
            call cont_value()

        ! Solve the HH Problem
            call solve_vfunction()
        
        ! Policy Functions - NOTE: COMMENT OUT WHEN DOING POLICY VS. COMPOSITION
            call policy()
        
        !*********************************************************** 
        
        ! **** THIS IS TO DECOMPOSE POLICY EFFECT ****
        ! **** FOR HH SIZE AND INCOME CHANNEL  ****
        ! **** before running: adjust section 'initialize' for respective counterfactual 
                 
        ! Policy function for risky share from counterfactual 
            !call policy1()  
           
        ! decomposition  --> return back to baseline 
            !call decomp1() 
            
        ! solve again continuation values during retirement
            !call cont_value()
            
        ! Solve again the HH Problem 
            !call solve_vfunction()
            
        ! All other policy function from baseline
            !call policy2() 
          
        !***********************************************************     
            
        ! ******  COMP. AND POL EFFECT - AGGREGATE ******** 
            
        ! Assign women the male policy function for risky share
            !call policy_gend()  
            
        ! All other policy functions as before
            !call policy2() 
            
        !***********************************************************     
         
        ! compute human capital returns: ONLY NEEDED WHEN COMPUTING HC RETURNS
	         ! call hc_return()
	       
	    !*********************************************************** 
         
        ! Simulate the economy
            call simulations() 
                        
        ! Export Results during working age
             call export()
            
        ! Export Results during retirement 
            !call export_ret()
            
         
contains


subroutine initialize

        ! deterministic income process by marital status and by gender
        ! note: 5000 is normalization factor in model, ihs transition of income
        
        ! constant 
        yfs = ( exp(9.6060d0+0.1375d0-0.7109d0) * exp(-log(2d0)) )/5000d0   
        yms = ( exp(9.6060d0+0.1375d0)* exp(-log(2d0)) ) /5000d0
        
        yfc = ( exp(7.5465d0+1.7196d0-3.9134d0) * exp(-log(2d0)) )/5000d0  
        ymc = ( exp(7.5465d0+1.7196d0)* exp(-log(2d0)) )/5000d0
        
        ! flat taxes (as in Kaplan & Huggett (2016))
        yfs = (1d0-tau)*yfs
        yms = (1d0-tau)*yms 
        
        yfc = (1d0-tau)*yfc
        ymc = (1d0-tau)*ymc
        
        ! age component   
        eff1fs = 0.0592d0
        eff2fs = -0.000728d0
        
        eff1ms = 0.0592d0
        eff2ms = -0.000728d0

        eff1fc = 0.0798d0
        eff2fc = -0.000884d0

        eff1mc = 0.0798d0
        eff2mc = -0.000884d0

		do j=1,JJ
        efffs(j)=exp((eff1fs+0.0110d0)*(j+minage) + eff2fs*(j+minage)**2d0)  
        effms(j)=exp(eff1ms*(j+minage) + eff2ms*(j+minage)**2d0)
        
        efffc(j)=exp((eff1fc+0.0523d0)*(j+minage) + eff2fc*(j+minage)**2d0)
        effmc(j)=exp(eff1mc*(j+minage) + eff2mc*(j+minage)**2d0)
        end do

        ! education component
        thetams(1)=exp(0.0d0)
        thetams(2)=exp(0.4480d0)
        thetafs(1)=exp(0.0d0)
        thetafs(2)=exp(0.4480d0+0.1160d0)
        
        thetamc(1)=exp(0.0d0)
        thetamc(2)=exp(0.5528d0)
        thetafc(1)=exp(0.0d0)
        thetafc(2)=exp(0.5528d0-0.3021d0)
        
        ! distribution of education; independent of marital status
        dist_thetaf(1)=0.44d0
        dist_thetaf(2)=0.56d0
        dist_thetam(1)=0.41d0
        dist_thetam(2)=0.59d0
        
        
       !! DECOMPOSITION EXERCISES: NEED TO COMMMENT OUT DEP. WHICH ONE TO RUN 
       
       ! decomposition 1
        
        !yfs     = yms
        !efffs   = effms
        !thetafs = thetams 
       
       ! decomposition 1b
       
        !yfs = yms
       
       ! decomposition 1c
       
        !efffs   = effms
       
       ! decomposition 1d
       
       ! thetafs = thetams
      
	   ! decomposition 1e (same level, but different slope)
	   ! check matlab code for explanation of individual steps
	   
	   ! efffs   = effms
	   ! yfs     = yfs*1.6679d0
	   
	   ! decomposition 1f (same slope, but different level)
	   ! check matlab code for explanation of individual steps

!	   do j=1,JJ
!	   add_low(j)  = ((yfs*efffs(j)*thetafs(1))+0.5954d0)/(yfs*efffs(j)*thetafs(1))
!	   add_high(j) = ((yfs*efffs(j)*thetafs(2))+0.5954d0)/(yfs*efffs(j)*thetafs(2))
!	   end do
	   
!	   ! adjustment of income process
!	   do j=1,JJ
!	   efffs(j) = efffs(j)*add_low(j)
!	   end do
!	   thetafs(2) = thetafs(2) *(add_high(10)/add_low(10))
	   
      ! decomposition 7
        
      ! dist_thetaf = dist_thetam

      ! stochstic part of income process by marital status and by gender
		
		! single women
        rhofs   	= 0.9348d0 
        varfs_boom  = 0.0856d0  
        varfs_rec   = 0.0856d0 
        mufs_rec    = -0.09d0  
        mufs_boom   = (-p_rec*mufs_rec)/(1d0-p_rec)
        
        rhoefs  = 0.0d0 
        varefs  = 0.1547d0 
        muefs   = 0.0d0
        
        ! married women 
        rhofc   	= 0.9273d0 
        varfc_boom  = 0.1583d0 
        varfc_rec   = 0.1583d0 
        mufc_rec    = -0.09d0 
        mufc_boom   = (-p_rec*mufc_rec)/(1d0-p_rec)
        
        rhoefc  = 0.0d0 
        varefc  = 0.2851d0 
        muefc   = 0.0d0
        
        ! single men 
        rhoms  	   = 0.9522d0 
        varms_boom = 0.0899d0  
        varms_rec  = 0.0899d0 
        mums_rec   = -0.09d0 
        mums_boom  = (-p_rec*mums_rec)/(1d0-p_rec)
        
        rhoems = 0.0d0 
        varems = 0.1689d0 
        muems  = 0.0d0
        
        ! married men 
        rhomc  	   = 0.9353d0 
        varmc_boom = 0.0809d0  
        varmc_rec  = 0.0809d0 
        mumc_rec   = -0.09d0  
        mumc_boom  = (-p_rec*mumc_rec)/(1d0-p_rec)
        
        rhoemc = 0.0d0 
        varemc = 0.1412d0 
        muemc  = 0.0d0
        
       ! decomposition 2
        !rhofs = rhoms
        !varfs_boom = varms_boom
        !varfs_rec = varms_rec
        !rhoefs = rhoems
        !varefs = varems
        
        
        ! discretize income shocks, persistent part
        
        call rouwenhorst(rhofs,sqrt(varfs_boom),mufs_boom,nz,zfs_boom,pizfs_boom)
        call rouwenhorst(rhofs,sqrt(varfs_rec),mufs_rec,nz,zfs_rec,pizfs_rec)
        
        call rouwenhorst(rhoms,sqrt(varms_boom),mums_boom,nz,zms_boom,pizms_boom) 
        call rouwenhorst(rhoms,sqrt(varms_rec),mums_rec,nz,zms_rec,pizms_rec) 
        
        ! transitory part (independent of aggregate state)
        
        call normal_discrete_1(efs, piefs(1,:), 0d0, varefs)
        call normal_discrete_1(ems, piems(1,:), 0d0, varems)
        
        do i=2,nz
        piefs(i,:) = piefs(1,:)
        piems(i,:) = piems(1,:)
        end do
        
        ! combining both shocks (depending on aggregate state)

		! single women: boom 
        do i=1,nz
        afs_boom((i-1)*nz+1:i*nz)=zfs_boom(i)+efs
        end do
        afs_boom=exp(afs_boom)
        
        ! single women: recession 
        do i=1,nz
        afs_rec((i-1)*nz+1:i*nz)=zfs_rec(i)+efs
        end do
        afs_rec=exp(afs_rec)
        
		! single men: boom
        do i=1,nz
        ams_boom((i-1)*nz+1:i*nz)=zms_boom(i)+ems
        end do
        ams_boom=exp(ams_boom)
        
        ! single men: recession
        do i=1,nz
        ams_rec((i-1)*nz+1:i*nz)=zms_rec(i)+ems
        end do
        ams_rec=exp(ams_rec)
            
        ! combined transtion matrix

		! single women: boom 
        do i=1,nz
        do j=1,nz
        pi1fs_boom(nz*(i-1)+1:nz*i,nz*(j-1)+1:nz*j)=pizfs_boom(i,j)
        end do
        end do
        
        ! single women: recession
        do i=1,nz
        do j=1,nz
        pi1fs_rec(nz*(i-1)+1:nz*i,nz*(j-1)+1:nz*j)=pizfs_rec(i,j)
        end do
        end do

		! transitory part (indpendent of aggregate state)
        do i=1,nz
        do j=1,nz
        pi2fs((nz*(i-1))+1:nz*i,(nz*(j-1))+1:nz*j)=piefs
        end do
        end do

		! single women: combined
        pifs_boom =pi1fs_boom*pi2fs
        pifs_rec  =pi1fs_rec*pi2fs
 
        
		! single men: boom 
        do i=1,nz
        do j=1,nz
        pi1ms_boom(nz*(i-1)+1:nz*i,nz*(j-1)+1:nz*j)=pizms_boom(i,j)
        end do
        end do
        
        ! single men: recession 
        do i=1,nz
        do j=1,nz
        pi1ms_rec(nz*(i-1)+1:nz*i,nz*(j-1)+1:nz*j)=pizms_rec(i,j)
        end do
        end do

		! transitory part (indpendent of aggregate state)
        do i=1,nz
        do j=1,nz
        pi2ms((nz*(i-1))+1:nz*i,(nz*(j-1))+1:nz*j)=piems
        end do
        end do

		! single men: combined
        pims_boom =pi1ms_boom*pi2ms
        pims_rec =pi1ms_rec*pi2ms
        
        
        ! matrix with income realizations in both states, needed for simulation
        afs_comb(:,1) = afs_boom
        afs_comb(:,2) = afs_rec
        
        ams_comb(:,1) = ams_boom
        ams_comb(:,2) = ams_rec
        
        
        ! income process of couples
        ! persistent part (depends on aggregate state)
        
        call rouwenhorst(rhofc,sqrt(varfc_boom),mufc_boom,nz,zfc_boom,pizfc_boom)
        call rouwenhorst(rhofc,sqrt(varfc_rec),mufc_rec,nz,zfc_rec,pizfc_rec)
        
        call rouwenhorst(rhomc,sqrt(varmc_boom),mumc_boom,nz,zmc_boom,pizmc_boom) 
        call rouwenhorst(rhomc,sqrt(varmc_rec),mumc_rec,nz,zmc_rec,pizmc_rec)

        ! transitory (independent of aggregate state but correlated across spouses)
        corr = 0.3d0  				! ROBUSTNESS: INCREASE TO 0.9d0
        sigmac = (/varefc,varemc/)
        mucoup = (/0d0, 0d0/)
        
        call normal_discrete_2(nn, coupgrid, Pcoup(1,:), mucoup, sigmac, corr)
        do i=2,nz*nz
        Pcoup(i,:) = Pcoup(1,:)
        end do

    
        ! grid of possible states within couple
        ! allocate persitent shock
    
        ! female component: boom 
        do j=1,nz
        do i=1,nzz*nz
        coupgrid_all_boom((j-1)*nzz*nz+i,1) =zfc_boom(j)
        end do
        end do
        
        ! female component: recession 
        do j=1,nz
        do i=1,nzz*nz
        coupgrid_all_rec((j-1)*nzz*nz+i,1) =zfc_rec(j)
        end do
        end do
        
        ! male component: boom 
        do k=1,nz
        do j=1,nz
        do i=1,nzz
        coupgrid_all_boom((k-1)*nzz*nz+(j-1)*nzz+i,2) =zmc_boom(j)
        end do
        end do
        end do
        
        ! male component: recession
        do k=1,nz
        do j=1,nz
        do i=1,nzz
        coupgrid_all_rec((k-1)*nzz*nz+(j-1)*nzz+i,2) =zmc_rec(j)
        end do
        end do
        end do
        
        ! allocate transitory shock
        do j= 1,nzz
        do i = 1,nzz
        coupgrid_all_boom((j-1)*nzz+i,:) = coupgrid_all_boom((j-1)*nzz+i,:)+ coupgrid(i,:)
        coupgrid_all_rec((j-1)*nzz+i,:)  =  coupgrid_all_rec((j-1)*nzz+i,:)+ coupgrid(i,:)
        end do
        end do 
   
        coupgrid_all_boom = exp(coupgrid_all_boom)
        coupgrid_all_rec  = exp(coupgrid_all_rec)
        
        
        ! one transition matrix within couple 
        ! transitory matrix (independent of aggregate state)
        do i=1,nzz
        do j=1,nzz
        pi1c((i-1)*nzz+1:nzz*i,(j-1)*nzz+1:nzz*j) = Pcoup
        end do 
        end do
        
        ! persistent part -- female, boom 
        do i=1,nz
        do j=1,nz
        pi2c_boom(nzz*nz*(i-1)+1:nzz*nz*i,nzz*nz*(j-1)+1:nzz*nz*j)=pizfc_boom(i,j)
        end do
        end do
        
        ! persistent part -- female, recession 
        do i=1,nz
        do j=1,nz
        pi2c_rec(nzz*nz*(i-1)+1:nzz*nz*i,nzz*nz*(j-1)+1:nzz*nz*j)=pizfc_rec(i,j)
        end do
        end do
        
        ! persistent part -- male, boom 
        do z=1,nz
        do p=1,nz
        do i=1,nz
        do j=1,nz
        pi3c_boom((z-1)*nzz*nz+nzz*(i-1)+1:(z-1)*nzz*nz+nzz*i,(p-1)*nzz*nz+nzz*(j-1)+1:(p-1)*nzz*nz+nzz*j)=pizmc_boom(i,j)
        end do
        end do
        end do
        end do
        
        ! persistent part -- male, recession 
        do z=1,nz
        do p=1,nz
        do i=1,nz
        do j=1,nz
        pi3c_rec((z-1)*nzz*nz+nzz*(i-1)+1:(z-1)*nzz*nz+nzz*i,(p-1)*nzz*nz+nzz*(j-1)+1:(p-1)*nzz*nz+nzz*j)=pizmc_rec(i,j)
        end do
        end do
        end do
        end do
       
        ! combining into one transtion matrix (depending on aggregate state)
        pic_boom = pi1c*pi2c_boom*pi3c_boom
        pic_rec  = pi1c*pi2c_rec*pi3c_rec
        
        
        ! return process risky asset - depends on aggregate state
        ! risky asset return during booms
        rhor      = 0.0d0
        varr 	  = 0.1758d0**2d0
        mur_boom  = 0.115d0 

        call rouwenhorst(rhor,sqrt(varr),mur_boom,nrr,ar_boom, pi_risk_boom)
        pir_boom = pi_risk_boom(1,:)
	
        ! risky asset return during recession (different mean, but same variance as in boom)
        mur_rec = -0.245d0 

        call rouwenhorst(rhor,sqrt(varr),mur_rec,nrr,ar_rec, pi_risk_rec)
        pir_rec = pi_risk_rec(1,:)
        
        ard_boom = 1d0+ar_boom
        ard_rec  = 1d0+ar_rec
  
        ! matrix that includes both states, needed for simulation
        ard_comb(:,1) = ard_boom
        ard_comb(:,2) = ard_rec
        
        ! grid for portfolio share (0 to 1) - equidistance
        risk_min = 0.0d0
        risk_max = 1.0d0
        call grid_Cons_Equi(rgrid,risk_min, risk_max)


        ! log spaced asset grid
        cmin = 0.0002d0  ! consumption floor
        cmax = 300.0d0
        c_lb=log(cmin+1.0d0)
        c_ub=log(cmax+1.0d0)

        call grid_Cons_Equi(cgrid,c_lb, c_ub)
        cgrid=exp(cgrid)-1

         
        ! equivalence scales 
           
        ! single men 
         num(1) = 1.115374d0
         num(2) = 1.116943d0
         num(3) = 1.125883d0
         num(4) = 1.156389d0
         num(5) = 1.179254d0
         num(6) = 1.198037d0
         num(7) = 1.177471d0
         num(8) = 1.156408d0
         num(9) = 1.155015d0
         num(10) = 1.149245d0
         num(11) = 1.197885d0
         num(12) = 1.182135d0
         num(13) = 1.181519d0
         num(14) = 1.143434d0
         num(15) = 1.179296d0
         num(16) = 1.159371d0
         num(17) = 1.185376d0
         num(18) = 1.15855d0
         num(19) = 1.184711d0
         num(20) = 1.204086d0
         num(21) = 1.201191d0
         num(22) = 1.175148d0
         num(23) = 1.162991d0
         num(24) = 1.164822d0
         num(25) = 1.140748d0
         num(26) = 1.145292d0
         num(27) = 1.127858d0
         num(28) = 1.117631d0
         num(29) = 1.202131d0
         num(30) = 1.129535d0
         num(31) = 1.175184d0
         num(32) = 1.095339d0
         num(33) = 1.142437d0
         num(34) = 1.077657d0
         num(35) = 1.116551d0
         num(36) = 1.069398d0

         ! single women 
         nuf(1) = 1.427366d0
         nuf(2) = 1.522981d0
         nuf(3) = 1.565953d0
         nuf(4) = 1.651393d0
         nuf(5) = 1.614807d0
         nuf(6) = 1.674927d0
         nuf(7) = 1.647139d0
         nuf(8) = 1.695039d0
         nuf(9) = 1.710022d0
         nuf(10) = 1.682682d0
         nuf(11) = 1.698723d0
         nuf(12) = 1.622424d0
         nuf(13) = 1.606513d0
         nuf(14) = 1.563697d0
         nuf(15) = 1.588745d0
         nuf(16) = 1.541049d0
         nuf(17) = 1.494798d0
         nuf(18) = 1.485049d0
         nuf(19) = 1.437968d0
         nuf(20) = 1.44561d0
         nuf(21) = 1.442471d0
         nuf(22) = 1.391488d0
         nuf(23) = 1.33112d0
         nuf(24) = 1.338498d0
         nuf(25) = 1.273979d0
         nuf(26) = 1.216328d0
         nuf(27) = 1.2448d0
         nuf(28) = 1.193775d0
         nuf(29) = 1.222786d0
         nuf(30) = 1.202131d0
         nuf(31) = 1.191707d0
         nuf(32) = 1.150265d0
         nuf(33) = 1.178141d0
         nuf(34) = 1.117523d0
         nuf(35) = 1.124148d0
         nuf(36) = 1.14483d0

        ! couples
         nuc(1) = 2.432881d0
         nuc(2) = 2.435853d0
         nuc(3) = 2.509238d0
         nuc(4) = 2.520984d0
         nuc(5) = 2.571975d0
         nuc(6) = 2.593857d0
         nuc(7) = 2.639374d0
         nuc(8) = 2.653389d0
         nuc(9) = 2.652963d0
         nuc(10) = 2.673717d0
         nuc(11) = 2.675789d0
         nuc(12) = 2.686544d0
         nuc(13) = 2.66684d0
         nuc(14) = 2.651449d0
         nuc(15) = 2.64894d0
         nuc(16) = 2.614885d0
         nuc(17) = 2.582108d0
         nuc(18) = 2.536446d0
         nuc(19) = 2.482672d0
         nuc(20) = 2.399675d0
         nuc(21) = 2.349969d0
         nuc(22) = 2.286824d0
         nuc(23) = 2.251347d0
         nuc(24) = 2.219593d0
         nuc(25) = 2.137313d0
         nuc(26) = 2.089818d0
         nuc(27) = 2.038475d0
         nuc(28) = 2.009359d0
         nuc(29) = 1.962437d0
         nuc(30) = 1.976248d0
         nuc(31) = 1.926663d0
         nuc(32) = 1.904638d0
         nuc(33) = 1.89657d0
         nuc(34) = 1.88797d0
         nuc(35) = 1.866706d0
         nuc(36) = 1.869405d0
         

         ! decomposition 3
         ! nuf = num

         ! transition matrix conditional on meeting a partner (seperately for men & women) 
         
         ! partner for men 
         trans_part_m(1,:) = (/0.08d0, 0.08d0, 0.08d0, 0.13d0, 0.13d0, 0.13d0, 0.05d0, 0.05d0, 0.05d0, &
                               0.02d0, 0.02d0, 0.02d0, 0.04d0, 0.04d0, 0.04d0, 0.02d0, 0.01d0, 0.01d0/)
         trans_part_m(2,:) = trans_part_m(1,:)    
         trans_part_m(3,:) = trans_part_m(1,:)
         
                               
         trans_part_m(4,:) = (/0.05d0, 0.05d0, 0.05d0, 0.11d0, 0.11d0, 0.11d0, 0.07d0, 0.07d0, 0.07d0, &
                               0.01d0, 0.01d0, 0.01d0, 0.05d0, 0.05d0, 0.05d0, 0.04d0, 0.04d0, 0.05d0/)                            
         trans_part_m(5,:) = trans_part_m(4,:)                            
         trans_part_m(6,:) = trans_part_m(4,:)
                               
                               
         trans_part_m(7,:) = (/0.05d0, 0.05d0, 0.05d0, 0.07d0, 0.07d0, 0.07d0, 0.12d0, 0.12d0, 0.12d0, &
                               0.01d0, 0.01d0, 0.01d0, 0.04d0, 0.04d0, 0.04d0, 0.05d0, 0.05d0, 0.03d0/)                            
         trans_part_m(8,:) = trans_part_m(7,:)                            
         trans_part_m(9,:) = trans_part_m(7,:)
                               
        
         trans_part_m(10,:) = (/0.03d0, 0.03d0, 0.03d0, 0.05d0, 0.05d0, 0.05d0, 0.02d0, 0.02d0, 0.02d0, &
                                0.05d0, 0.05d0, 0.05d0, 0.11d0, 0.11d0, 0.11d0, 0.08d0, 0.08d0, 0.06d0/)
         trans_part_m(11,:) = trans_part_m(10,:)                           
         trans_part_m(12,:) = trans_part_m(10,:)
                
                                
         trans_part_m(13,:) = (/0.04d0, 0.04d0, 0.04d0, 0.07d0, 0.07d0, 0.07d0, 0.00d0, 0.00d0, 0.00d0, &
                                0.09d0, 0.09d0, 0.09d0, 0.14d0, 0.14d0, 0.12d0, 0.00d0, 0.00d0, 0.00d0/)                             
         trans_part_m(14,:) = trans_part_m(13,:)                        
         trans_part_m(15,:) = trans_part_m(13,:)
                                
         
         trans_part_m(16,:) = (/0.01d0, 0.01d0, 0.01d0, 0.02d0, 0.02d0, 0.02d0, 0.07d0, 0.07d0, 0.07d0, &
                                0.05d0, 0.05d0, 0.05d0, 0.07d0, 0.07d0, 0.07d0, 0.11d0, 0.11d0, 0.12d0/)                                   
         trans_part_m(17,:) = trans_part_m(16,:)
         trans_part_m(18,:) = trans_part_m(16,:)
         
         
         ! partner for women      
         trans_part_f(1,:) = (/0.06d0, 0.06d0, 0.06d0, 0.09d0, 0.09d0, 0.09d0, 0.06d0, 0.06d0, 0.06d0, &
                               0.04d0, 0.04d0, 0.04d0, 0.06d0, 0.06d0, 0.06d0, 0.02d0, 0.02d0, 0.03d0/)                         
         trans_part_f(2,:) = trans_part_f(1,:)                     
         trans_part_f(3,:) = trans_part_f(1,:)
         
                                
         trans_part_f(4,:) = (/0.06d0, 0.06d0, 0.06d0, 0.10d0, 0.10d0, 0.10d0, 0.05d0, 0.05d0, 0.05d0, &
                               0.04d0, 0.04d0, 0.04d0, 0.06d0, 0.06d0, 0.06d0, 0.02d0, 0.02d0, 0.03d0/)                      
         trans_part_f(5,:) =  trans_part_f(4,:)                      
         trans_part_f(6,:) =  trans_part_f(4,:)

                               
         trans_part_f(7,:) = (/0.03d0, 0.03d0, 0.03d0, 0.09d0, 0.09d0, 0.09d0, 0.11d0, 0.11d0, 0.11d0, &
                               0.02d0, 0.02d0, 0.02d0, 0.00d0, 0.00d0, 0.00d0, 0.090d0, 0.08d0, 0.08d0/)
         trans_part_f(8,:) =  trans_part_f(7,:)                     
         trans_part_f(9,:) =  trans_part_f(7,:)
                               
                               
         trans_part_f(10,:) = (/0.01d0, 0.01d0, 0.01d0, 0.02d0, 0.02d0, 0.02d0, 0.01d0, 0.01d0, 0.01d0, &
                                0.06d0, 0.06d0, 0.06d0, 0.12d0, 0.12d0, 0.12d0, 0.11d0, 0.11d0, 0.12d0/)                               
         trans_part_f(11,:) = trans_part_f(10,:)
         trans_part_f(12,:) = trans_part_f(10,:)
                               
                               
         trans_part_f(13,:) = (/0.02d0, 0.02d0, 0.02d0, 0.04d0, 0.04d0, 0.04d0, 0.03d0, 0.03d0, 0.03d0, &
                                0.07d0, 0.07d0, 0.07d0, 0.10d0, 0.10d0, 0.10d0, 0.07d0, 0.07d0, 0.08d0/)                         
         trans_part_f(14,:) = trans_part_f(13,:)         
         trans_part_f(15,:) = trans_part_f(13,:)
                               
                               
         trans_part_f(16,:) = (/0.01d0, 0.01d0, 0.01d0, 0.05d0, 0.05d0, 0.05d0, 0.05d0, 0.05d0, 0.05d0, &
                                0.07d0, 0.070d0, 0.07d0, 0.00d0, 0.00d0, 0.00d0, 0.15d0, 0.15d0, 0.16d0/)                                
         trans_part_f(17,:) = trans_part_f(16,:)
         trans_part_f(18,:) = trans_part_f(16,:)
             
    
         ! marriage and divorce probabilities
        
         ! low education
         do l=1,JJ
         mu_m(l,1,:) = exp(0.1202d0-0.0588d0*(l+minage))/ &
                  (1d0+exp(0.1202d0-0.0588d0*(l+minage)))

         mu_f(l,1,:) = exp(0.1202d0-0.0588d0*(l+minage)-0.3203d0)/ &
                  (1d0+exp(0.1202d0-0.0588d0*(l+minage)-0.3203d0))
         end do
         
         ! high education
         do l=1,JJ
         mu_m(l,2,:) = exp(0.1202d0-0.0588d0*(l+minage)+ 0.0896d0)/ &
                  (1d0+exp(0.1202d0-0.0588d0*(l+minage)+ 0.0896d0))
         
         mu_f(l,2,:) =  exp(0.1202d0-0.0588d0*(l+minage)-0.3203d0+ 0.0896d0)/ &
                   (1d0+exp(0.1202d0-0.0588d0*(l+minage)-0.3203d0+ 0.0896d0))
         end do
         
         
         
         ! education: low-low
         
            do l=1,JJ
                    do zi=1,nzz
                        do zp=1,nzz
         lambda(l,1,1,zi,zp)=exp(-1.5753d0-0.0442d0*(l+minage))/ &
                        (1d0+exp(-1.5753d0-0.0442d0*(l+minage)))
                        end do
                    end do
            end do
            
        ! education: low-high
        
            do l=1,JJ
                    do zi=1,nzz
                        do zp=1,nzz
         lambda(l,1,2,zi,zp)=exp(-1.5753d0-0.0442d0*(l+minage)-0.4061d0)/ &
                        (1d0+exp(-1.5753d0-0.0442d0*(l+minage)-0.4061d0))
                    end do
                end do
            end do
            
        ! education: high-high
        
            do l=1,JJ
                    do zi=1,nzz
                        do zp=1,nzz
        lambda(l,2,2,zi,zp)=exp(-1.5753d0-0.0442d0*(l+minage)-0.0826d0-0.4061d0)/ &
                       (1d0+exp(-1.5753d0-0.0442d0*(l+minage)-0.0826d0-0.4061d0))
                    end do
                end do
            end do
            
        ! education: high-low
        
            do l=1,JJ
                    do zi=1,nzz
                        do zp=1,nzz
         lambda(l,2,1,zi,zp)=exp(-1.5753d0-0.0442d0*(l+minage)-0.0826d0)/ &
                        (1d0+exp(-1.5753d0-0.0442d0*(l+minage)-0.0826d0))
                    end do
                end do
            end do


        ! decomposition 4
        !  mu_f = mu_m


        ! average net worth by gender, age and education (for marriage market)
		! men, low edu
         mwealth(1,1,:) = (/ 6.37d0, 6.66d0, 6.74d0, 7.21d0, 11.01d0, 7.67d0, &
                            16.25d0, 5.74d0, 10.04d0, 6.99d0, 21.20d0, 12.66d0, &
                            9.65d0, 14.12d0, 8.78d0, 13.43d0, 13.25d0, 17.10d0, &
                            29.43d0, 14.83d0, 25.45d0, 15.63d0, 30.26d0, 16.88d0, &
                            32.55d0, 18.32d0, 35.42d0, 17.84d0, 19.03d0, 38.24d0, &
                            33.23d0, 50.86d0, 16.71d0, 26.44d0, 15.07d0, 25.05d0, 25.05d0 /) 
         
         ! men, high edu
         mwealth(1,2,:) = (/ 13.32d0, 16.85d0, 20.13d0, 10.05d0, 13.85d0, 33.08d0, &
                             19.53d0, 29.69d0, 36.75d0, 41.34d0, 26.18d0, 42.46d0, &
                             38.29d0, 37.38d0, 38.03d0, 68.41d0, 40.34d0, 62.51d0, &
                             55.55d0, 61.45d0, 82.37d0, 110.72d0, 121.48d0, 75.83d0, &
                             137.02d0, 79.10d0, 88.43d0, 110.07d0, 78.65d0, 127.51d0, &
                             101.52d0, 105.96d0, 112.82d0, 168.16d0, 111.42d0, 69.65d0, 69.65d0 /)
         
         ! women, low edu
         mwealth(2,1,:) = (/ 2.07d0, 2.30d0, 1.27d0, 2.89d0, 3.99d0, 3.78d0, &
                             3.42d0, 3.69d0, 6.21d0, 3.90d0, 5.97d0, 8.59d0, &
                             9.86d0, 4.52d0, 8.76d0, 12.95d0, 8.91d0, 8.23d0, &
                             8.15d0, 14.92d0, 6.08d0, 12.85d0, 11.01d0, 14.36d0, &
                             17.29d0, 15.40d0, 11.96d0, 16.67d0, 20.18d0, 18.33d0, &
                             21.13d0, 16.17d0, 31.86d0, 29.26d0, 20.84d0, 19.52d0, 19.52d0 /)  
         
         ! women, high edu
         mwealth(2,2,:) = (/ 6.02d0, 4.43d0, 11.54d0, 7.91d0, 6.13d0, 15.34d0, &
                             14.0d0, 8.67d0, 9.32d0, 9.59d0, 20.38d0, 16.50d0, &
                             19.19d0, 26.74d0, 49.74d0, 22.0d0, 38.14d0, 31.02d0, &
                             42.05d0, 39.16d0, 52.33d0, 46.25d0, 62.12d0, 47.87d0, &
                             49.21d0, 52.93d0, 79.08d0, 45.57d0, 64.37d0, 67.35d0, &
                             72.92d0, 51.43d0, 74.63d0, 78.19d0, 67.23d0, 68.14d0, 68.14d0 /) 
             
                             
         ! decomposition 5
         ! trans_part_f = trans_part_m       ! decomposition 5a
         ! mwealth(1,:,:) = mwealth(2,:,:)	 ! decomposition 5b
         

        ! export parameters of the income process
         
        ! deterministic part
         open(unit=1,file='results/yfs.txt', &
         status='replace',action='write')
         open(unit=2,file='results/yms.txt', &
         status='replace',action='write')
         open(unit=3,file='results/yfc.txt', &
         status='replace',action='write')
         open(unit=4,file='results/ymc.txt', &
         status='replace',action='write')
         open(unit=5,file='results/efffs.txt', &
         status='replace',action='write')
         open(unit=6,file='results/effms.txt', &
         status='replace',action='write')
         open(unit=7,file='results/efffc.txt', &
         status='replace',action='write')
         open(unit=8,file='results/effmc.txt', &
         status='replace',action='write')
         open(unit=9,file='results/thetafs.txt', &
         status='replace',action='write')
         open(unit=10,file='results/thetams.txt', &
         status='replace',action='write')
         open(unit=11,file='results/thetafc.txt', &
         status='replace',action='write')
         open(unit=12,file='results/thetamc.txt', &
         status='replace',action='write')
         
        ! stochastic part
         open(unit=13,file='results/afs_boom.txt', &
         status='replace',action='write')
		 open(unit=14,file='results/afs_rec.txt', &
         status='replace',action='write')
         open(unit=15,file='results/ams_boom.txt', &
         status='replace',action='write')
         open(unit=16,file='results/ams_rec.txt', &
         status='replace',action='write')
         open(unit=17,file='results/coupgrid_all_boom.txt', &
         status='replace',action='write')
         open(unit=18,file='results/coupgrid_all_rec.txt', &
         status='replace',action='write')
         open(unit=19,file='results/pic_boom.txt', &
         status='replace',action='write')         
         open(unit=20,file='results/pic_rec.txt', &
         status='replace',action='write')
         open(unit=21,file='results/pifs_boom.txt', &
         status='replace',action='write')         
         open(unit=22,file='results/pifs_rec.txt', &
         status='replace',action='write')
         open(unit=23,file='results/pims_boom.txt', &
         status='replace',action='write')         
         open(unit=24,file='results/pims_rec.txt', &
         status='replace',action='write')
         open(unit=25,file='results/ard_boom.txt', &
         status='replace',action='write')         
         open(unit=26,file='results/ard_rec.txt', &
         status='replace',action='write')
         
         
        write(1,*) yfs
        write(2,*) yms
        write(3,*) yfc
        write(4,*) ymc
        write(5,*) efffs
        write(6,*) effms
        write(7,*) efffc
        write(8,*) effmc
        write(9,*)  thetafs
        write(10,*) thetams
        write(11,*) thetafc
        write(12,*) thetamc 
        write(13,*) afs_boom
        write(14,*) afs_rec
        write(15,*) ams_boom
        write(16,*) ams_rec
        write(17,*) coupgrid_all_boom
        write(18,*) coupgrid_all_rec
        write(19,*) pic_boom
        write(20,*) pic_rec
        write(21,*) pifs_boom
        write(22,*) pifs_rec
        write(23,*) pims_boom
        write(24,*) pims_rec
        write(25,*) ard_boom
        write(26,*) ard_rec
         
		do i=1,26
			close(unit=i)
		end do
       
  write(*,*)'initialization done'
         
end subroutine

subroutine cont_value

        ! see file retirement.f90

           call initialize_ret()
           call solve_vfret()
           call policy_ret()
           call VF_retirement()
           call export_ret()
	       !call hcval_ret()  ! comment in when computing human capital return


end subroutine

subroutine solve_vfunction

    
    do l=JJ,1,-1            ! age loop (has to be first loop)
        
	    if (l==JJ) then
			ll= l
	    else
			ll = l+1
        end if
        
       ! ******** parallelize ********
        
       !$OMP PARALLEL PRIVATE(maxbndm,maxbndc,maxbndf,minbndc,minbndf,minbndm,tp,zp,p,output_m,output_f,output_c,t,z)
              
       !$OMP DO 
       
       do i=1,ncc

        
        ! start loop over remaining state variables
        
        do t=1,naa
        
               do z=1,nzz

					do p=1,nrisk 
			
        ! ******** single male ********
        
        ! specify bounds for optimization (HH chooses cprime, state variable is asset cash-in-hand)
        minbndm = cmin
        if(p == 1) then
        maxbndm = max(0.0d0,min(cgrid(i),cmax))       
        else
        maxbndm = max(0.0d0,min(cgrid(i)-SMFs(l),cmax))
        end if
        
        
        ! optimization
        output_m = &
        golden_sing(minbndm,maxbndm,rgrid(p),gammam,yms,thetams(t),effms(ll),&
					betam,cgrid(i),l,t,z,VF_m(:,:,t,ll),VF_cm(:,:,:,t,:,ll),&
                    pims_boom,pims_rec,mu_m(l,t,z),trans_part_m,num(l),cont_m(:,z,t),SMFs(l), &
                    mwealth(2,:,l),ams_boom,ams_rec,ymc,yfc,effmc(ll),efffc(ll),thetamc(t),thetafc,1)
                    

        cash_opt_m(p,i,z,t,l)    = output_m(1)  ! this is cprime
        Vopt_m(p,i,z,t,l)        = output_m(2)
        copt_m(p,i,z,t,l)        = output_m(3)
        ctopt_sing_m(p,i,z,t,l)  = output_m(4)
        ctopt_coup_m(p,i,z,t,l)  = output_m(5)
        
        
        ! ******** single female ********
       
        ! specify bounds for optimization
        minbndf = cmin
        if(p == 1) then 
        maxbndf = max(0.0d0,min(cgrid(i),cmax))
        else
        maxbndf = max(0.0d0,min(cgrid(i)-SMFs(l),cmax))
        end if
        
        ! optimization
        output_f = &
        golden_sing(minbndf,maxbndf,rgrid(p),gammaf,yfs,thetafs(t),efffs(ll),&
					betaf,cgrid(i),l,t,z,VF_f(:,:,t,ll),VF_cf(:,:,:,t,:,ll),&
                    pifs_boom,pifs_rec,mu_f(l,t,z),trans_part_f,nuf(l),cont_f(:,z,t),SMFs(l), &
                    mwealth(1,:,l),afs_boom,afs_rec,yfc,ymc,efffc(ll),effmc(ll),thetafc(t),thetamc,2)


        cash_opt_f(p,i,z,t,l)    = output_f(1) ! this is cprime
        Vopt_f(p,i,z,t,l)        = output_f(2)
        copt_f(p,i,z,t,l)        = output_f(3)
        ctopt_sing_f(p,i,z,t,l)  = output_f(4)
        ctopt_coup_f(p,i,z,t,l)  = output_f(5)
        
        
        
       ! ******** couples  ********
    
              do tp=1, naa
              do zp=1,nzz
			  
			  
        ! specify bounds for optimization
        minbndc = cmin
        if(p == 1) then
        maxbndc = max(0.0d0,min(cgrid(i),cmax))                      
        else
        maxbndc = max(0.0d0,min(cgrid(i)-SMFc(l),cmax))
        end if
        
        ! optimization
        output_c = &
        golden_coup(minbndc,maxbndc,rgrid(p),gammac,gammac,yfc,ymc, &
                    thetafc(t),thetamc(tp),efffc(ll),effmc(ll),betac,betac, &
                    cgrid(i),l,z,zp,VF_f(:,:,t,ll), VF_m(:,:,tp,ll),VF_c(:,:,:,t,tp,ll), &
                    VF_cf(:,:,:,t,tp,ll),VF_cm(:,:,:,tp,t,ll),pifs_boom, pifs_rec, &
                    pims_boom,pims_rec,pic_boom,pic_rec,lambda(l,t,tp,z,zp),nuc(l), &
                    cont_c(:,zp,t,tp),cont_ci(:,zp,t,tp),cont_cp(:,zp,t,tp),SMFc(l),yfs,yms,&
                    efffs(ll),effms(ll),thetafs(t),thetams(tp),afs_boom,afs_rec,ams_boom,ams_rec)
        
        ! output    
        cash_opt_c(p,i,z,zp,t,tp,l)  = output_c(1) ! this is cprime
        Vopt_c(p,i,z,zp,t,tp,l)      = output_c(2)
        Vcoup_f(p,i,z,zp,t,tp,l)     = output_c(3)
        Vcoup_m(p,i,zp,z,tp,t,l)     = output_c(4)
        copt_c(l,t,tp,i,z,zp,p)      = output_c(5)
        
  
                end do    
                end do
                
                
        end do    ! end of nrisk loop      
        end do    ! end of nzz loop
        end do    ! end of naa loop
        end do    ! end of ncc loop          
                        
            !$OMP END DO
            !$OMP END PARALLEL 
      
      
          
        ! Find optimal share of risky assets 
        
        do i = 1,ncc
            do t=1,naa
                do z= 1,nzz
         
        ! singles men  
        inx2 = 0
        inx  = 0       
                    
        VF_m(i,z,t,l)   = maxval(Vopt_m(:,i,z,t,l))
        inx2            = maxloc(Vopt_m(:,i,z,t,l))
        inx             = maxval(inx2,1)
        

        ! policy function for consumption
        cpol_m(i,z,t,l)    = copt_m(inx,i,z,t,l)
        
        ! continuation value couple/ single
        ct_sing_m(i,z,t,l) = ctopt_sing_m(inx,i,z,t,l)
        ct_coup_m(i,z,t,l) = ctopt_coup_m(inx,i,z,t,l)
          
        ! singles women 
        inx2 = 0
        inx  = 0 
                    
        VF_f(i,z,t,l)   = maxval(Vopt_f(:,i,z,t,l))
        inx2            = maxloc(Vopt_f(:,i,z,t,l))
        inx             = maxval(inx2,1)
        
        ! policy function for  consumption
        cpol_f(i,z,t,l)    = copt_f(inx,i,z,t,l)
        
        ! continuation value couple/ single
        ct_sing_f(i,z,t,l) = ctopt_sing_f(inx,i,z,t,l)
        ct_coup_f(i,z,t,l) = ctopt_coup_f(inx,i,z,t,l)
        
                end do
            end do
        end do
        
        ! couples
        do i = 1,ncc
            do t=1,naa
                do z= 1,nzz
                
                  do tp=1, naa
                  do zp=1,nzz
                  
        inx2 = 0
        inx  = 0 

        VF_c(i,z,zp,t,tp,l)  = maxval(Vopt_c(:,i,z,zp,t,tp,l))
        
        inx2 = maxloc(Vopt_c(:,i,z,zp,t,tp,l))
        inx  = maxval(inx2,1)
        
        ! Policy function for individual in couple      
        VF_cf(i,z,zp,t,tp,l)      = Vcoup_f(inx,i,z,zp,t,tp,l)
        VF_cm(i,zp,z,tp,t,l)      = Vcoup_m(inx,i,zp,z,tp,t,l)
        
        ! policy function for consumption
        cpol_c(l,t,tp,i,z,zp)     = copt_c(l,t,tp,i,z,zp,inx)
        
                    end do
                    end do
                
                end do
            end do
        end do
         
    
       write(*,'(a,i3,a)')'Age: ',l+minage,' DONE!'
        
           
        end do         ! end of age loop 
        
end subroutine

subroutine policy


     ! ASSET RELATED POLICY FUNCTIONS


     ! single women  
        do l=1,JJ
            do t=1,naa
                do i=1,ncc
                    do z=1,nzz
                    
        inx2 = 0
        inx  = 0 
                    
        inx2 = maxloc(Vopt_f(:,i,z,t,l))
        inx  = maxval(inx2,1)
        
        cprimepol_f(i,z,t,l) = cash_opt_f(inx,i,z,t,l)
        sharepol_f(i,z,t,l)  = rgrid(inx)
        riskypol_f(i,z,t,l)  = sharepol_f(i,z,t,l)*cprimepol_f(i,z,t,l) 
        savepol_f(i,z,t,l)   = (1d0-sharepol_f(i,z,t,l))*cprimepol_f(i,z,t,l)
        
        
                    end do
                end do
            end do
        end do


        ! single men
        
        do l=1,JJ
            do t=1,naa
                do i=1,ncc
                    do z=1,nzz
                    
        inx2 = 0
        inx  = 0 
                    
        inx2 = maxloc(Vopt_m(:,i,z,t,l))
        inx  = maxval(inx2,1)
                    
        cprimepol_m(i,z,t,l) = cash_opt_m(inx,i,z,t,l)
        sharepol_m(i,z,t,l)  = rgrid(inx)
        riskypol_m(i,z,t,l)  = sharepol_m(i,z,t,l)*cprimepol_m(i,z,t,l) 
        savepol_m(i,z,t,l)   = (1d0-sharepol_m(i,z,t,l))*cprimepol_m(i,z,t,l)

                    end do
                end do
            end do
        end do
        

        ! couple
        
        do l=1,JJ
            do ti=1,naa
              do tp=1,naa
                do i=1,ncc
                    do zi=1,nzz
                      do zp=1,nzz
                      
        inx2 = 0
        inx  = 0 
                      
        inx2 = maxloc(Vopt_c(:,i,zi,zp,ti,tp,l))
        inx  = maxval(inx2,1)
        
                    
        cprimepol_c(l,ti,tp,i,zi,zp) = cash_opt_c(inx,i,zi,zp,ti,tp,l)
        sharepol_c(l,ti,tp,i,zi,zp)  = rgrid(inx)
        riskypol_c(l,ti,tp,i,zi,zp)  = sharepol_c(l,ti,tp,i,zi,zp)*cprimepol_c(l,ti,tp,i,zi,zp) 
        savepol_c(l,ti,tp,i,zi,zp)   = (1d0-sharepol_c(l,ti,tp,i,zi,zp))*cprimepol_c(l,ti,tp,i,zi,zp)

                      end do
                    end do
                end do
              end do
            end do
        end do
        
        
    write(*,*) 'policy functions done'
    

end subroutine

subroutine hc_return

!single women, low education [assume never get married]
call comp_hc_val(savepol_f, riskypol_f, yfs, efffs, thetafs(1), afs_boom, afs_rec, pifs_boom, pifs_rec,cpol_f,1,nuf, &
				distret_f_low, cpolret_f, pen_f(:,1), m_yf,m_efff,m_thetaf(1), psif, dist_f_low, hc_val_f_low) 


! export results
!open(unit=1,file='results/dist_f_low.txt', &
!status='replace', action='write')
!open(unit=2,file='results/hc_val_f_low.txt', &
!status='replace', action='write')

!write(1,*)dist_f_low
!write(2,*)hc_val_f_low

!do i=1,2
!close(unit=i)
!end do 

!single women, high education [assume never get married]
call comp_hc_val(savepol_f, riskypol_f, yfs, efffs, thetafs(2), afs_boom, afs_rec, pifs_boom, pifs_rec, cpol_f,2,nuf, &
				distret_f_high, cpolret_f, pen_f(:,2), m_yf,m_efff,m_thetaf(2),psif,dist_f_high, hc_val_f_high) 

!! export results
!open(unit=1,file='results/dist_f_high.txt', &
!status='replace', action='write')
!open(unit=2,file='results/hc_val_f_high.txt', &
!status='replace', action='write')

!write(1,*)dist_f_high
!write(2,*)hc_val_f_high

!do i=1,2
!close(unit=i)
!end do 

!single men, low education [assume never get married] 
call comp_hc_val(savepol_m, riskypol_m, yms, effms, thetams(1), ams_boom, ams_rec, pims_boom, pims_rec, cpol_m,1,num, &
				distret_m_low, cpolret_m, pen_m(:,1), m_ym,m_effm,m_thetam(1), psim, dist_m_low, hc_val_m_low) 

! export results
!open(unit=1,file='results/dist_m_low.txt', &
!status='replace', action='write')
!open(unit=2,file='results/hc_val_m_low.txt', &
!status='replace', action='write')

!write(1,*)dist_m_low
!write(2,*)hc_val_m_low

!do i=2,2
!close(unit=i)
!end do 

!single men, high education [assume never get married] 
call comp_hc_val(savepol_m, riskypol_m, yms, effms, thetams(2), ams_boom, ams_rec, pims_boom, pims_rec, cpol_m,2,num, &
				distret_m_high, cpolret_m, pen_m(:,2), m_ym,m_effm,m_thetam(2), psim, dist_m_high, hc_val_m_high) 

!open(unit=1,file='results/dist_m_high.txt', &
!status='replace', action='write')
!open(unit=2,file='results/hc_val_m_high.txt', &
!status='replace', action='write')

!write(1,*)dist_m_high
!write(2,*)hc_val_m_high

!do i=2,2
!close(unit=i)
!end do 


end subroutine

subroutine comp_hc_val(safepol, riskypol,ys, eff, thetas, as_boom, as_rec,pi_boom, pi_rec, &			!computing hc values
						cpol, ind, nu, distret, cpolret,pen,m_y,m_eff,m_theta, surv,dist, hc_val)

real*8, intent(in)  :: safepol(ncc,nzz,naa,JJ), riskypol(ncc,nzz,naa,JJ), cpol(ncc,nzz,naa,JJ), nu(JJ)

real*8, intent(in)  :: ys, eff(JJ+1), thetas, as_boom(nzz), as_rec(nzz), pi_boom(nzz,nzz), pi_rec(nzz,nzz) 

real*8, intent(in)  :: distret(JR,nzz*ncc,nzz*ncc), cpolret(ncc,nzz,naa,JR)

real*8, intent(in)  :: pen(nzz), m_y, m_eff(JR), m_theta, surv(JR)

integer, intent(in) :: ind

real*8, intent(out) :: dist(JJ,nzz*ncc,nzz*ncc), hc_val(JJ,ncc,nzz)


integer :: gridl_boom, gridl_rec, gridh_boom, gridh_rec

integer :: l, i, zz, zpr, apr, rpr, kret, k

real*8, dimension(JJ,ncc,nzz,nzz,nrr) :: apr_boom, apr_rec
        
real*8  :: diff_boom, diff_rec, auxi,  prob2(nzz),  prob3(nzz,nzz), prob2_ret(nzz)

real*8  :: hlp_boom, hlp_rec , sdf(JJ,ncc,nzz,JJ,ncc,nzz), sdf_ret(JJ,ncc,nzz,JR,ncc,nzz), surv_dis

real*8  :: hc_val_ret(JJ,ncc,nzz), auxi_ret, prob3_ret(nzz,nzz), help
        
real*8, dimension(nzz*ncc,nzz*ncc) :: prob, prob_ret

! compute transition matrix across state space except for terminal period
dist = 0d0

do l = 1, JJ-1
	do i = 1,ncc
		do zz = 1,nzz
		
	do zpr = 1,nzz    ! next period values (prime)
	do rpr = 1,nrr
	
!determine aprime given shock outcomes, depending on aggregrate state
apr_boom(l,i,zz,zpr,rpr) = &
(1d0+r)*safepol(i,zz,ind,l) + ys*eff(l+1)*thetas*as_boom(zpr) + (1d0+ar_boom(rpr))*riskypol(i,zz,ind,l)

apr_rec(l,i,zz,zpr,rpr) = &
(1d0+r)*safepol(i,zz,ind,l) + ys*eff(l+1)*thetas*as_rec(zpr)  + (1d0+ar_rec(rpr))*riskypol(i,zz,ind,l)

! translate into point on asset grid
gridl_boom = max(cntarray(cgrid,ncc,apr_boom(l,i,zz,zpr,rpr)),1)   
gridh_boom = min(gridl_boom+1,ncc)

gridl_rec = max(cntarray(cgrid,ncc,apr_rec(l,i,zz,zpr,rpr)),1)   
gridh_rec = min(gridl_rec+1,ncc)

if (gridl_boom .NE. gridh_boom) then 
diff_boom = max((apr_boom(l,i,zz,zpr,rpr)- cgrid(gridl_boom))/(cgrid(gridh_boom)-cgrid(gridl_boom)),0d0)
elseif (gridl_boom == gridh_boom) then
diff_boom = 0d0 
end if 

if (gridl_rec .NE. gridh_rec) then 
diff_rec  = max((apr_rec(l,i,zz,zpr,rpr) - cgrid(gridl_rec)) /(cgrid(gridh_rec)-cgrid(gridl_rec)),0d0)
elseif (gridl_rec == gridh_rec) then
diff_rec = 0d0 
end if 

! always allocate to higher grid point (to avoid houshols ending up in lowest grid point with positive probability)
diff_boom = 1d0
diff_rec = 1d0

! fill matrix of probabilities
dist(l,(zz-1)*ncc+i,(zpr-1)*ncc+gridh_boom) = dist(l,(zz-1)*ncc+i,(zpr-1)*ncc+gridh_boom) + diff_boom*pi_boom(zz,zpr)*(1d0-p_rec)*pir_boom(rpr)
dist(l,(zz-1)*ncc+i,(zpr-1)*ncc+gridl_boom) = dist(l,(zz-1)*ncc+i,(zpr-1)*ncc+gridl_boom) + (1d0-diff_boom)*pi_boom(zz,zpr)*(1d0-p_rec)*pir_boom(rpr)

dist(l,(zz-1)*ncc+i,(zpr-1)*ncc+gridh_rec) = dist(l,(zz-1)*ncc+i,(zpr-1)*ncc+gridh_rec) + diff_rec*pi_rec(zz,zpr)*p_rec*pir_rec(rpr)
dist(l,(zz-1)*ncc+i,(zpr-1)*ncc+gridl_rec) = dist(l,(zz-1)*ncc+i,(zpr-1)*ncc+gridl_rec) + (1d0-diff_rec)*pi_rec(zz,zpr)*p_rec*pir_rec(rpr)

	 end do
	 end do

			end do
		end do
	end do
	
! compute transition matrix for terminal period
	do i = 1,ncc
		do zz = 1,nzz
			do rpr = 1,nrr
	
apr_boom(JJ,i,zz,zz,rpr) = &
(1d0+r)*safepol(i,zz,ind,JJ) + pen(zz)-(m_y*m_eff(1)*m_theta) + (1d0+ar_boom(rpr))*riskypol(i,zz,ind,JJ)

apr_rec(JJ,i,zz,zz,rpr) = &
(1d0+r)*safepol(i,zz,ind,JJ) + pen(zz)-(m_y*m_eff(1)*m_theta) + (1d0+ar_rec(rpr))*riskypol(i,zz,ind,JJ)


! translate into point on asset grid
gridl_boom = max(cntarray(cgrid,ncc,apr_boom(JJ,i,zz,zz,rpr)),1)   
gridh_boom = min(gridl_boom+1,ncc)

gridl_rec = max(cntarray(cgrid,ncc,apr_rec(JJ,i,zz,zz,rpr)),1)   
gridh_rec = min(gridl_rec+1,ncc)

if (gridl_boom .NE. gridh_boom) then 
diff_boom = max((apr_boom(JJ,i,zz,zz,rpr)- cgrid(gridl_boom))/(cgrid(gridh_boom)-cgrid(gridl_boom)),0d0)
elseif (gridl_boom == gridh_boom) then
diff_boom = 0d0 
end if 

if (gridl_rec .NE. gridh_rec) then 
diff_rec  = max((apr_rec(JJ,i,zz,zz,rpr) - cgrid(gridl_rec)) /(cgrid(gridh_rec)-cgrid(gridl_rec)),0d0)
elseif (gridl_rec == gridh_rec) then
diff_rec = 0d0 
end if 

! always allocate to higher grid point (to avoid houshols ending up in lowest grid point with positive probability)
diff_boom = 1d0
diff_rec = 1d0

! fill matrix of probabilities
dist(JJ,(zz-1)*ncc+i,(zz-1)*ncc+gridh_boom) = dist(JJ,(zz-1)*ncc+i,(zz-1)*ncc+gridh_boom) + diff_boom*(1d0-p_rec)*pir_boom(rpr)
dist(JJ,(zz-1)*ncc+i,(zz-1)*ncc+gridl_boom) = dist(JJ,(zz-1)*ncc+i,(zz-1)*ncc+gridl_boom) + (1d0-diff_boom)*(1d0-p_rec)*pir_boom(rpr)

dist(JJ,(zz-1)*ncc+i,(zz-1)*ncc+gridh_rec) = dist(JJ,(zz-1)*ncc+i,(zz-1)*ncc+gridh_rec) + diff_rec*p_rec*pir_rec(rpr)
dist(JJ,(zz-1)*ncc+i,(zz-1)*ncc+gridl_rec) = dist(JJ,(zz-1)*ncc+i,(zz-1)*ncc+gridl_rec) + (1d0-diff_rec)*p_rec*pir_rec(rpr)

			
			end do
		end do
	end do
	
		
! compute the stochastic discount factor except for terminal period
sdf = 0d0
sdf_ret = 0d0
do l = 1, JJ-1 

prob = 0d0

! remaining periods during working life
	do k = l+1,JJ

prob = dist(l,:,:)

if (k > l+1) then
do j =l+2,k
prob = MATMUL(prob,dist(j-1,:,:))
end do
end if

do i = 1,ncc
	do zz = 1, nzz
		do apr = 1, ncc
		do zpr = 1, nzz	
		
if (prob((zz-1)*ncc+i,(zpr-1)*ncc+apr)==0) then
sdf(l,i,zz,k,apr,zpr) = 0d0
else
sdf(l,i,zz,k,apr,zpr) = (((nu(k)**gammaf)*cpol(apr,zpr,ind,k)**(-gammaf))/((nu(l)**gammaf)*cpol(i,zz,ind,l)**(-gammaf))) &
*(betaf**k)*(prob((zz-1)*ncc+i,(zpr-1)*ncc+apr))
end if

		end do
		end do
	 end do
end do

	 end do ! end of k loop 
	 
	 
! remaining periods during retirement
do kret = 1, JR

! generate prob_ret for first period of retirement 
prob_ret = dist(l,:,:)
do j =l+1,JJ
prob_ret = MATMUL(prob_ret,dist(j,:,:))
end do

! expand for remaining retirement periods
if (kret > 1) then
do j =2,kret
prob_ret = MATMUL(prob_ret,distret(j-1,:,:))
end do
end if

! uncertain survival during retirment
surv_dis = surv(1)
if (kret > 1) then
do j =2,kret
surv_dis = surv_dis * surv(j)
end do
end if


do i = 1,ncc
	do zz = 1, nzz
		do apr = 1, ncc
		do zpr = 1, nzz
			
if (prob_ret((zz-1)*ncc+i,(zpr-1)*ncc+apr)==0) then
sdf_ret(l,i,zz,kret,apr,zpr) = 0d0
else
sdf_ret(l,i,zz,kret,apr,zpr) = ((cpolret(apr,zpr,ind,kret)**(-gammaf))/((nu(l)**gammaf)*cpol(i,zz,ind,l)**(-gammaf)))* &
prob_ret((zz-1)*ncc+i,(zpr-1)*ncc+apr)*(betaf**(JJ-l+kret))*surv_dis 
end if

		end do
		end do
	 end do
end do 

	end do ! end of kret loop 

end do ! end of age loop 


! SDF in terminal working period (JJ)
do kret = 1, JR

prob_ret = dist(JJ,:,:)
if (kret > 1) then
do j =2,kret
prob_ret = MATMUL(prob_ret,distret(j-1,:,:))
end do
end if

! uncertain survival during retirment
surv_dis = surv(1)
if (kret > 1) then
do j =2,kret
surv_dis = surv_dis * surv(j)
end do
end if



do i = 1,ncc
	do zz = 1, nzz
		do apr = 1, ncc
	
if (prob_ret((zz-1)*ncc+i,(zz-1)*ncc+apr)==0) then
sdf(JJ,i,zz,kret,apr,zz) = 0d0
else
sdf(JJ,i,zz,kret,apr,zz) = ((cpolret(apr,zz,ind,kret)**(-gammaf))/((nu(JJ)**gammaf)*cpol(i,zz,ind,JJ)**(-gammaf)))* &
prob_ret((zz-1)*ncc+i,(zz-1)*ncc+apr)*(betaf**kret)*surv_dis 
end if

		end do
	 end do
end do 

end do ! end of kret loop (terminal period)



! compute human capital values during working life 
hc_val = 0d0
do l = 1,JJ-1
	do zz = 1, nzz
		do k = l+1,JJ

auxi = 0d0
prob2 = 0d0
		
! compute expected earnings
if (k == l+1) then
hlp_boom = 0d0
hlp_rec = 0d0
do j = 1,nzz
hlp_boom = hlp_boom + pi_boom(zz,j)*as_boom(j)
hlp_rec =  hlp_rec  + pi_rec(zz,j)*as_rec(j)
end do
auxi = (1d0-p_rec)*hlp_boom + p_rec*hlp_rec

elseif (k > l+1) then
prob2 = (1d0-p_rec)*pi_boom(zz,:)+p_rec*pi_rec(zz,:)
prob3 = (1d0-p_rec)*pi_boom+p_rec*pi_rec

if (k > l+2) then
do j =1,k-(l+2) 
prob2 = MATMUL(prob2,prob3)
end do
end if

hlp_boom = 0d0
hlp_rec = 0d0
do j = 1,nzz
hlp_boom = hlp_boom + prob2(j)*as_boom(j)
hlp_rec =  hlp_rec  + prob2(j)*as_rec(j)
end do
auxi = (1d0-p_rec)*hlp_boom + p_rec*hlp_rec
end if

do i = 1, ncc
help = sum(sdf(l,i,zz,k,:,:))*ys*eff(k)*thetas*auxi
hc_val(l,i,zz) = hc_val(l,i,zz) + help
end do


		end do
	end do
end do

! compute human capital values during retirement
hc_val_ret = 0d0
do l = 1,JJ-1
	do zz = 1, nzz

! expected prod state in period JJ (to determine pensions)
prob2_ret = (1d0-p_rec)*pi_boom(zz,:)+p_rec*pi_rec(zz,:)
prob3_ret = (1d0-p_rec)*pi_boom+p_rec*pi_rec

if (l<JJ-1) then
do j =l+1,JJ 
prob2_ret = MATMUL(prob2_ret,prob3_ret)
end do
end if

auxi_ret = 0d0
do j = 1,nzz
auxi_ret = auxi_ret + prob2_ret(j)*pen(j)
end do 


do kret = 1,JR
	do i = 1, ncc
help=sum(sdf_ret(l,i,zz,kret,:,:))*(auxi_ret-(m_y*m_eff(kret)*m_theta))
hc_val_ret(l,i,zz) = hc_val_ret(l,i,zz) + help	
	end do			
end do	
		
	end do
end do

! compute human capital value in terminal period 
do zz = 1, nzz
	do kret = 1,JR
		do i = 1, ncc
help = sum(sdf(JJ,i,zz,kret,:,:))*(pen(zz)-(m_y*m_eff(kret)*m_theta))		
hc_val(JJ,i,zz) = hc_val(JJ,i,zz) + help		
		end do			
	end do	
end do


! combine working life and retirement (do not need to do it for terminal working period)
do l = 1,JJ-1
	do zz = 1, nzz
		do i = 1, ncc
		
hc_val(l,i,zz) = hc_val(l,i,zz) + hc_val_ret(l,i,zz)

		end do
	end do			
end do	



end subroutine comp_hc_val

subroutine simulations

        ! see file simulate_economy.f90

            call shock_path()
            call simulate()
            call cohort_vars()

end subroutine

subroutine export

        ! store results in txt files (singles)

        open(unit=1,file='results/vff.txt', &
        status='replace', action='write')
        open(unit=2,file='results/vfm.txt', &
        status='replace', action='write')
        open(unit=3,file='results/share_f.txt', &
        status='replace',action='write')
        open(unit=4,file='results/share_m.txt', &
        status='replace',action='write')
        open(unit=5,file='results/cprime_f.txt', &
        status='replace',action='write')
        open(unit=6,file='results/cprime_m.txt', &
        status='replace',action='write')
        open(unit=7,file='results/risky_f.txt', &
        status='replace',action='write')
        open(unit=8,file='results/risky_m.txt', &
        status='replace',action='write')
        open(unit=9,file='results/safe_f.txt', &
        status='replace',action='write')
        open(unit=10,file='results/safe_m.txt', &
        status='replace',action='write')
        open(unit=11,file='results/consum_f.txt', &
        status='replace',action='write')
        open(unit=12,file='results/consum_m.txt', &
        status='replace',action='write')
        open(unit=13,file='results/ct_sing_f.txt', &
        status='replace',action='write')
        open(unit=14,file='results/ct_coup_f.txt', &
        status='replace',action='write')
        open(unit=15,file='results/ct_sing_m.txt', &
        status='replace',action='write')
        open(unit=16,file='results/ct_coup_m.txt', &
        status='replace',action='write')
        
        write(1,*)VF_f
        write(2,*)VF_m
        write(3,*)sharepol_f
        write(4,*)sharepol_m
        write(5,*)cprimepol_f
        write(6,*)cprimepol_m
        write(7,*)riskypol_f
        write(8,*)riskypol_m
        write(9,*)savepol_f
        write(10,*)savepol_m
        write(11,*)cpol_f
        write(12,*)cpol_m
        write(13,*)ct_sing_f
        write(14,*)ct_coup_f
        write(15,*)ct_sing_m
        write(16,*)ct_coup_m
        
        ! couples
        open(unit=17,file='results/vfc.txt', &
        status='replace',action='write')
        open(unit=18,file='results/share_c.txt', &
        status='replace',action='write')
        open(unit=19,file='results/cprime_c.txt', &
        status='replace',action='write')
        open(unit=20,file='results/risky_c.txt', &
        status='replace',action='write')
        open(unit=21,file='results/safe_c.txt', &
        status='replace',action='write')
        open(unit=22,file='results/consum_c.txt', &
        status='replace',action='write')
        
        write(17,*)VF_c
        write(18,*)sharepol_c
        write(19,*)cprimepol_c
        write(20,*)riskypol_c
        write(21,*)savepol_c
        write(22,*)cpol_c
     
        
        ! cohort variables
        open(unit=23,file='results/cohort/cashcoh_sf.txt', &
        status='replace',action='write')
        open(unit=24,file='results/cohort/consumcoh_sf.txt', &
        status='replace',action='write')
        open(unit=25,file='results/cohort/sharecoh_sf.txt', &
        status='replace',action='write')
        open(unit=26,file='results/cohort/riskycoh_sf.txt', &
        status='replace',action='write')
        open(unit=27,file='results/cohort/safecoh_sf.txt', &
        status='replace',action='write')
        open(unit=28,file='results/cohort/SMPcoh_sf.txt', &
        status='replace',action='write')
        open(unit=29,file='results/cohort/ucsharecoh_sf.txt', &
        status='replace',action='write')
        
        
        open(unit=30,file='results/cohort/cashcoh_sm.txt', &
        status='replace',action='write')
        open(unit=31,file='results/cohort/consumcoh_sm.txt', &
        status='replace',action='write')
        open(unit=32,file='results/cohort/sharecoh_sm.txt', &
        status='replace',action='write')
        open(unit=33,file='results/cohort/riskycoh_sm.txt', &
        status='replace',action='write')
        open(unit=34,file='results/cohort/safecoh_sm.txt', &
        status='replace',action='write')
        open(unit=35,file='results/cohort/SMPcoh_sm.txt', &
        status='replace',action='write')
        open(unit=36,file='results/cohort/ucsharecoh_sm.txt', &
        status='replace',action='write')
        
        
        open(unit=37,file='results/cohort/cashcoh_c.txt', &
        status='replace',action='write')
        open(unit=38,file='results/cohort/consumcoh_c.txt', &
        status='replace',action='write')
        open(unit=39,file='results/cohort/sharecoh_c.txt', &
        status='replace',action='write')
        open(unit=40,file='results/cohort/riskycoh_c.txt', &
        status='replace',action='write')
        open(unit=41,file='results/cohort/safecoh_c.txt', &
        status='replace',action='write')
        open(unit=42,file='results/cohort/SMPcoh_c.txt', &
        status='replace',action='write')
        open(unit=43,file='results/cohort/ucsharecoh_c.txt', &
        status='replace',action='write')


        open(unit=44,file='results/VF_cf.txt', &
        status='replace',action='write')
        open(unit=45,file='results/VF_cm.txt', &
        status='replace',action='write')
        open(unit=46,file='results/Vopt_f.txt', &
        status='replace',action='write')
        
        open(unit=47,file='results/cohort/hcval_coh_sf.txt', &
        status='replace',action='write')
        open(unit=48,file='results/cohort/hcret_coh_sf.txt', &
        status='replace',action='write')
        open(unit=49,file='results/cohort/hcval_coh_sm.txt', &
        status='replace',action='write')
        open(unit=50,file='results/cohort/hcret_coh_sm.txt', &
        status='replace',action='write')
        
        
        write(23,*)cash_coh_sf
        write(24,*)c_coh_sf
        write(25,*)share_coh_sf
        write(26,*)risky_coh_sf
        write(27,*)safe_coh_sf
        write(28,*)SMP_coh_sf
        write(29,*)ucshare_coh_sf
        
        write(30,*)cash_coh_sm
        write(31,*)c_coh_sm
        write(32,*)share_coh_sm
        write(33,*)risky_coh_sm
        write(34,*)safe_coh_sm
        write(35,*)SMP_coh_sm
        write(36,*)ucshare_coh_sm
        
        write(37,*)cash_coh_c
        write(38,*)c_coh_c
        write(39,*)share_coh_c
        write(40,*)risky_coh_c
        write(41,*)safe_coh_c
        write(42,*)SMP_coh_c
        write(43,*)ucshare_coh_c
        
        write(44,*)VF_cf
        write(45,*)VF_cm
        write(46,*)Vopt_f
        
        write(47,*)hcval_coh_sf
        write(48,*)hcret_coh_sf
        write(49,*)hcval_coh_sm
        write(50,*)hcret_coh_sm
        
        do i=1,50
        close(unit=i)
        end do 
        
        ! simulation variables 
        open(unit=51,file='results/simdata/chain_martm.txt', &
        status='replace',action='write')
        open(unit=52,file='results/simdata/chain_martf.txt', &
        status='replace',action='write')
        open(unit=53,file='results/simdata/chain_prodm.txt', &
        status='replace',action='write')
        open(unit=54,file='results/simdata/chain_prodf.txt', &
        status='replace',action='write')
        open(unit=55,file='results/simdata/chain_prodmm.txt', &
        status='replace',action='write')
        open(unit=56,file='results/simdata/chain_prodmf.txt', &
        status='replace',action='write')
        open(unit=57,file='results/simdata/assetgrid.txt', &
        status='replace',action='write')
        open(unit=58,file='results/simdata/grid_path_m.txt', &
        status='replace',action='write')
        open(unit=59,file='results/simdata/grid_path_f.txt', &
        status='replace',action='write')
        open(unit=60,file='results/simdata/grid_path_c.txt', &
        status='replace',action='write')
        open(unit=61,file='results/simdata/share_path_sm.txt', &
        status='replace',action='write')
        open(unit=62,file='results/simdata/share_path_sf.txt', &
        status='replace',action='write')
        open(unit=63,file='results/simdata/share_path_c.txt', &
        status='replace',action='write')
        open(unit=64,file='results/simdata/SMP_path_sm.txt', &
        status='replace',action='write')
        open(unit=65,file='results/simdata/SMP_path_sf.txt', &
        status='replace',action='write')
        open(unit=66,file='results/simdata/SMP_path_c.txt', &
        status='replace',action='write')
        open(unit=67,file='results/simdata/ams_boom.txt', &
        status='replace',action='write')
        open(unit=68,file='results/simdata/ams_rec.txt', &
        status='replace',action='write')
        open(unit=69,file='results/simdata/afs_boom.txt', &
        status='replace',action='write')
        open(unit=70,file='results/simdata/afs_rec.txt', &
        status='replace',action='write')
        open(unit=71,file='results/simdata/cprime_path_c.txt', &
        status='replace',action='write')
        open(unit=72,file='results/simdata/gridl_path_f.txt', &
        status='replace',action='write')
        open(unit=73,file='results/simdata/gridh_path_f.txt', &
        status='replace',action='write')
        open(unit=74,file='results/simdata/grid_path_f_t.txt', &
        status='replace',action='write')
        open(unit=75,file='results/simdata/cntsing_path_sf.txt', &
        status='replace',action='write')
        open(unit=76,file='results/simdata/cntcoup_path_sf.txt', &
        status='replace',action='write')
        open(unit=77,file='results/simdata/cntsing_path_sm.txt', &
        status='replace',action='write')
        open(unit=78,file='results/simdata/cntcoup_path_sm.txt', &
        status='replace',action='write')
        open(unit=79,file='results/simdata/hcval_path_sf.txt', &
        status='replace',action='write')
        open(unit=80,file='results/simdata/hcret_path_sf.txt', &
        status='replace',action='write')
        open(unit=81,file='results/simdata/hcval_path_sm.txt', &
        status='replace',action='write')
        open(unit=82,file='results/simdata/hcret_path_sm.txt', &
        status='replace',action='write')
        open(unit=83,file='results/simdata/c_path_sf.txt', &
        status='replace',action='write')
        open(unit=84,file='results/simdata/c_path_sm.txt', &
        status='replace',action='write')
        
        
        write(51,*)chain_martm
        write(52,*)chain_martf
        write(53,*)chain_prodm
        write(54,*)chain_prodf
        write(55,*)chain_prodmm
        write(56,*)chain_prodmf
        write(57,*)cgrid
        write(58,*)cash_path_sm
        write(59,*)cash_path_sf
        write(60,*)cash_path_c
        write(61,*)share_path_sm
        write(62,*)share_path_sf
        write(63,*)share_path_c
        write(64,*)SMP_path_sm
        write(65,*)SMP_path_sf
        write(66,*)SMP_path_c
        write(67,*)ams_boom
        write(68,*)ams_rec
        write(69,*)afs_boom
        write(70,*)afs_rec
        write(71,*)cprime_path_c
        write(72,*)gridl_path_f
        write(73,*)gridh_path_f
        write(74,*)grid_path_f
        write(75,*)cntsing_path_sf
        write(76,*)cntcoup_path_sf
        write(77,*)cntsing_path_sm
        write(78,*)cntcoup_path_sm
        write(79,*)hcval_path_sf
        write(80,*)hcret_path_sf
        write(81,*)hcval_path_sm
        write(82,*)hcret_path_sm
        write(83,*)c_path_sf
        write(84,*)c_path_sm


   do i=50,84
        close(unit=i)
        end do
       
    
end subroutine

subroutine decomp1

		! values from baseline

		!!!! DETERMINISTIC INCOME COMPONENT
	    ! constant 
         yfs = ( exp(9.6060d0+0.1375d0-0.7109d0) * exp(-log(2d0)) )/5000d0 
        
        ! flat taxes (as in Kaplan & Huggett (2016))
        yfs = (1d0-tau)*yfs
        
        ! age component   
        eff1fs = 0.0592d0
        eff2fs = -0.000728d0
		do j=1,JJ
        efffs(j)=exp((eff1fs+0.0110d0)*(j+minage) + eff2fs*(j+minage)**2d0)  
        end do

        ! education component
        thetafs(1)=exp(0.0d0)
        thetafs(2)=exp(0.4480d0+0.1160d0)

		!!!! HOUSEHOLD SIZES
        ! single women 
         nuf(1) = 1.427366d0
         nuf(2) = 1.522981d0
         nuf(3) = 1.565953d0
         nuf(4) = 1.651393d0
         nuf(5) = 1.614807d0
         nuf(6) = 1.674927d0
         nuf(7) = 1.647139d0
         nuf(8) = 1.695039d0
         nuf(9) = 1.710022d0
         nuf(10) = 1.682682d0
         nuf(11) = 1.698723d0
         nuf(12) = 1.622424d0
         nuf(13) = 1.606513d0
         nuf(14) = 1.563697d0
         nuf(15) = 1.588745d0
         nuf(16) = 1.541049d0
         nuf(17) = 1.494798d0
         nuf(18) = 1.485049d0
         nuf(19) = 1.437968d0
         nuf(20) = 1.44561d0
         nuf(21) = 1.442471d0
         nuf(22) = 1.391488d0
         nuf(23) = 1.33112d0
         nuf(24) = 1.338498d0
         nuf(25) = 1.273979d0
         nuf(26) = 1.216328d0
         nuf(27) = 1.2448d0
         nuf(28) = 1.193775d0
         nuf(29) = 1.222786d0
         nuf(30) = 1.202131d0
         nuf(31) = 1.191707d0
         nuf(32) = 1.150265d0
         nuf(33) = 1.178141d0
         nuf(34) = 1.117523d0
         nuf(35) = 1.124148d0
         nuf(36) = 1.14483d0
        
!!		! STOCHASTIC INCOME COMPONENT
		! single women
        rhofs   	= 0.9348d0 
        varfs_boom  = 0.0856d0  
        varfs_rec   = 0.0856d0 
        mufs_rec    = -0.09d0  
        mufs_boom   = (-p_rec*mufs_rec)/(1d0-p_rec)
        
        rhoefs  = 0.0d0 
        varefs  = 0.1547d0 
        muefs   = 0.0d0
        
       ! discretize income shocks
        call rouwenhorst(rhofs,sqrt(varfs_boom),mufs_boom,nz,zfs_boom,pizfs_boom)
        call rouwenhorst(rhofs,sqrt(varfs_rec),mufs_rec,nz,zfs_rec,pizfs_rec)
             
        ! transitory part (independent of aggregate state)
        call normal_discrete_1(efs, piefs(1,:), 0d0, varefs)
        
        do i=2,nz
        piefs(i,:) = piefs(1,:)
        end do
        
        ! combining both shocks (depending on aggregate state)

		! single women: boom 
        do i=1,nz
        afs_boom((i-1)*nz+1:i*nz)=zfs_boom(i)+efs
        end do
        afs_boom=exp(afs_boom)
        
        ! single women: recession 
        do i=1,nz
        afs_rec((i-1)*nz+1:i*nz)=zfs_rec(i)+efs
        end do
        afs_rec=exp(afs_rec)
      
        ! combined transtion matrix

		! single women: boom 
        do i=1,nz
        do j=1,nz
        pi1fs_boom(nz*(i-1)+1:nz*i,nz*(j-1)+1:nz*j)=pizfs_boom(i,j)
        end do
        end do
        
        ! single women: recession
        do i=1,nz
        do j=1,nz
        pi1fs_rec(nz*(i-1)+1:nz*i,nz*(j-1)+1:nz*j)=pizfs_rec(i,j)
        end do
        end do

		! transitory part (indpendent of aggregate state)
        do i=1,nz
        do j=1,nz
        pi2fs((nz*(i-1))+1:nz*i,(nz*(j-1))+1:nz*j)=piefs
        end do
        end do

		! single women: combined
        pifs_boom =pi1fs_boom*pi2fs
        pifs_rec  =pi1fs_rec*pi2fs       
        
        ! matrix with income realizations in both states, needed for simulation
        afs_comb(:,1) = afs_boom
        afs_comb(:,2) = afs_rec

end subroutine

subroutine policy1


     ! SINGLE WOMEN POLICY FOR RISKY SHARE FROM COUNTERFACTUAL SIMULATION

     ! single women    
          do l=1,JJ
            do t=1,naa
                do i=1,ncc
                    do z=1,nzz
                    
        inx2 = 0
        inx  = 0 
                    
        inx2 = maxloc(Vopt_f(:,i,z,t,l))
        inx  = maxval(inx2,1)
        
        sharepol_f(i,z,t,l)  = rgrid(inx)
        
                    end do
                end do
            end do
        end do
    
end subroutine

subroutine policy2

     ! ALL REMAINING POLICY FUNCTIONS (EXCEPT RISKY ASSET SINGLE WOMEN)
     
       ! single men  
        do l=1,JJ
            do t=1,naa
                do i=1,ncc
                    do z=1,nzz
                    
        inx2 = 0
        inx  = 0 
                    
        inx2 = maxloc(Vopt_m(:,i,z,t,l))
        inx  = maxval(inx2,1)
                    
        sharepol_m(i,z,t,l)  = rgrid(inx)

                    end do
                end do
            end do
        end do
        

        ! couple
        do l=1,JJ
            do ti=1,naa
              do tp=1,naa
                do i=1,ncc
                    do zi=1,nzz
                      do zp=1,nzz
                      
        inx2 = 0
        inx  = 0 
                      
        inx2 = maxloc(Vopt_c(:,i,zi,zp,ti,tp,l))
        inx  = maxval(inx2,1)
                    
        sharepol_c(l,ti,tp,i,zi,zp)  = rgrid(inx)

                      end do
                    end do
                end do
              end do
            end do
        end do
        
        
     
     
     ! ASSET RELATED POLICY FUNCTIONS

     ! single women 
        do l=1,JJ
            do t=1,naa
                do i=1,ncc
                    do z=1,nzz
                    
        inx2 = 0
        inx  = 0 
                    
        inx2 = maxloc(Vopt_f(:,i,z,t,l))
        inx  = maxval(inx2,1)
        
        cprimepol_f(i,z,t,l) = cash_opt_f(inx,i,z,t,l)
        riskypol_f(i,z,t,l)  = sharepol_f(i,z,t,l)*cprimepol_f(i,z,t,l) 
        savepol_f(i,z,t,l)   = (1d0-sharepol_f(i,z,t,l))*cprimepol_f(i,z,t,l)
        
        
                    end do
                end do
            end do
        end do


      ! single men 
        do l=1,JJ
            do t=1,naa
                do i=1,ncc
                    do z=1,nzz
                    
        inx2 = 0
        inx  = 0 
                    
        inx2 = maxloc(Vopt_m(:,i,z,t,l))
        inx  = maxval(inx2,1)
                    
        cprimepol_m(i,z,t,l) = cash_opt_m(inx,i,z,t,l)
        riskypol_m(i,z,t,l)  = sharepol_m(i,z,t,l)*cprimepol_m(i,z,t,l) 
        savepol_m(i,z,t,l)   = (1d0-sharepol_m(i,z,t,l))*cprimepol_m(i,z,t,l)

                    end do
                end do
            end do
        end do
        
        
        
	  ! couples
        do l=1,JJ
            do ti=1,naa
              do tp=1,naa
                do i=1,ncc
                    do zi=1,nzz
                      do zp=1,nzz
                      
        inx2 = 0
        inx  = 0 
                      
        inx2 = maxloc(Vopt_c(:,i,zi,zp,ti,tp,l))
        inx  = maxval(inx2,1)
        
        cprimepol_c(l,ti,tp,i,zi,zp) = cash_opt_c(inx,i,zi,zp,ti,tp,l)
        riskypol_c(l,ti,tp,i,zi,zp)  = sharepol_c(l,ti,tp,i,zi,zp)*cprimepol_c(l,ti,tp,i,zi,zp) 
        savepol_c(l,ti,tp,i,zi,zp)   = (1d0-sharepol_c(l,ti,tp,i,zi,zp))*cprimepol_c(l,ti,tp,i,zi,zp)

                      end do
                    end do
                end do
              end do
            end do
        end do
    
end subroutine

subroutine policy_gend

    ! ASSIGN SINGLE WOMEN THE MALE POLICIES FOR THE RISKY SHARE

	! policy function risky share for single men 
        
        do l=1,JJ
            do t=1,naa
                do i=1,ncc
                    do z=1,nzz
                    
        inx2 = 0
        inx  = 0 
                    
        inx2 = maxloc(Vopt_m(:,i,z,t,l))
        inx  = maxval(inx2,1)
                    
         sharepol_m(i,z,t,l)  = rgrid(inx)

                    end do
                end do
            end do
        end do
        
        
       ! assign women the same policy for risky share as men 
        
          do l=1,JJ
            do t=1,naa
                do i=1,ncc
                    do z=1,nzz
                    
        sharepol_f(i,z,t,l)  =  sharepol_m(i,z,t,l)
        
                    end do
                end do
            end do
        end do



end subroutine

end program main
